home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / linalg / test.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  50.6 KB  |  1,829 lines
  1. /* linalg/test.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman, Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* Author:  G. Jungman
  21.  */
  22. #include <config.h>
  23. #include <stdlib.h>
  24. #include <gsl/gsl_test.h>
  25. #include <gsl/gsl_math.h>
  26. #include <gsl/gsl_ieee_utils.h>
  27. #include <gsl/gsl_permute_vector.h>
  28. #include <gsl/gsl_blas.h>
  29. #include <gsl/gsl_complex_math.h>
  30. #include <gsl/gsl_linalg.h>
  31.  
  32. int check (double x, double actual, double eps);
  33. gsl_matrix * create_hilbert_matrix(size_t size);
  34. gsl_matrix * create_general_matrix(size_t size1, size_t size2);
  35. gsl_matrix * create_vandermonde_matrix(size_t size);
  36. gsl_matrix * create_moler_matrix(size_t size);
  37. gsl_matrix * create_row_matrix(size_t size1, size_t size2);
  38. gsl_matrix * create_2x2_matrix(double a11, double a12, double a21, double a22);
  39. int test_matmult(void);
  40. int test_matmult_mod(void);
  41. int test_LU_solve_dim(const gsl_matrix * m, const double * actual, double eps);
  42. int test_LU_solve(void);
  43. int test_LUc_solve_dim(const gsl_matrix_complex * m, const double * actual, double eps);
  44. int test_LUc_solve(void);
  45. int test_QR_solve_dim(const gsl_matrix * m, const double * actual, double eps);
  46. int test_QR_solve(void);
  47. int test_QR_lssolve_dim(const gsl_matrix * m, const double * actual, double eps);
  48. int test_QR_lssolve(void);
  49. int test_QR_decomp_dim(const gsl_matrix * m, double eps);
  50. int test_QR_decomp(void);
  51. int test_QRPT_solve_dim(const gsl_matrix * m, const double * actual, double eps);
  52. int test_QRPT_solve(void);
  53. int test_QRPT_decomp_dim(const gsl_matrix * m, double eps);
  54. int test_QRPT_decomp(void);
  55. int test_QR_update_dim(const gsl_matrix * m, double eps);
  56. int test_QR_update(void);
  57. int test_SV_solve_dim(const gsl_matrix * m, const double * actual, double eps);
  58. int test_SV_solve(void);
  59. int test_SV_decomp_dim(const gsl_matrix * m, double eps);
  60. int test_SV_decomp(void);
  61. int test_cholesky_solve_dim(const gsl_matrix * m, const double * actual, double eps);
  62. int test_cholesky_solve(void);
  63. int test_cholesky_decomp_dim(const gsl_matrix * m, double eps);
  64. int test_cholesky_decomp(void);
  65. int test_HH_solve_dim(const gsl_matrix * m, const double * actual, double eps);
  66. int test_HH_solve(void);
  67. int test_TD_solve_dim(size_t dim, double d, double od, const double * actual, double eps);
  68. int test_TD_solve(void);
  69. int test_bidiag_decomp_dim(const gsl_matrix * m, double eps);
  70. int test_bidiag_decomp(void);
  71.  
  72. int 
  73. check (double x, double actual, double eps)
  74. {
  75.   if (actual == 0)
  76.     return fabs(x) > eps;
  77.   else
  78.     return (fabs(x - actual)/fabs(actual)) > eps;
  79. }
  80.  
  81. gsl_matrix *
  82. create_hilbert_matrix(size_t size)
  83. {
  84.   size_t i, j;
  85.   gsl_matrix * m = gsl_matrix_alloc(size, size);
  86.   for(i=0; i<size; i++) {
  87.     for(j=0; j<size; j++) {
  88.       gsl_matrix_set(m, i, j, 1.0/(i+j+1.0));
  89.     }
  90.   }
  91.   return m;
  92. }
  93.  
  94. gsl_matrix *
  95. create_general_matrix(size_t size1, size_t size2)
  96. {
  97.   size_t i, j;
  98.   gsl_matrix * m = gsl_matrix_alloc(size1, size2);
  99.   for(i=0; i<size1; i++) {
  100.     for(j=0; j<size2; j++) {
  101.       gsl_matrix_set(m, i, j, 1.0/(i+j+1.0));
  102.     }
  103.   }
  104.   return m;
  105. }
  106.  
  107. gsl_matrix *
  108. create_singular_matrix(size_t size1, size_t size2)
  109. {
  110.   size_t i, j;
  111.   gsl_matrix * m = gsl_matrix_alloc(size1, size2);
  112.   for(i=0; i<size1; i++) {
  113.     for(j=0; j<size2; j++) {
  114.       gsl_matrix_set(m, i, j, 1.0/(i+j+1.0));
  115.     }
  116.   }
  117.  
  118.   /* zero the first column */
  119.   for(j = 0; j <size2; j++)
  120.     gsl_matrix_set(m,0,j,0.0);
  121.  
  122.   return m;
  123. }
  124.  
  125.  
  126. gsl_matrix *
  127. create_vandermonde_matrix(size_t size)
  128. {
  129.   size_t i, j;
  130.   gsl_matrix * m = gsl_matrix_alloc(size, size);
  131.   for(i=0; i<size; i++) {
  132.     for(j=0; j<size; j++) {
  133.       gsl_matrix_set(m, i, j, pow(i + 1.0, size - j - 1.0));
  134.     }
  135.   }
  136.   return m;
  137. }
  138.  
  139. gsl_matrix *
  140. create_moler_matrix(size_t size)
  141. {
  142.   size_t i, j;
  143.   gsl_matrix * m = gsl_matrix_alloc(size, size);
  144.   for(i=0; i<size; i++) {
  145.     for(j=0; j<size; j++) {
  146.       gsl_matrix_set(m, i, j, GSL_MIN(i+1,j+1)-2.0);
  147.     }
  148.   }
  149.   return m;
  150. }
  151.  
  152. gsl_matrix_complex *
  153. create_complex_matrix(size_t size)
  154. {
  155.   size_t i, j;
  156.   gsl_matrix_complex * m = gsl_matrix_complex_alloc(size, size);
  157.   for(i=0; i<size; i++) {
  158.     for(j=0; j<size; j++) {
  159.       gsl_complex z = gsl_complex_rect(1.0/(i+j+1.0), 1/(i*i+j*j+0.5));
  160.       gsl_matrix_complex_set(m, i, j, z);
  161.     }
  162.   }
  163.   return m;
  164. }
  165.  
  166. gsl_matrix *
  167. create_row_matrix(size_t size1, size_t size2)
  168. {
  169.   size_t i;
  170.   gsl_matrix * m = gsl_matrix_calloc(size1, size2);
  171.   for(i=0; i<size1; i++) {
  172.       gsl_matrix_set(m, i, 0, 1.0/(i+1.0));
  173.   }
  174.  
  175.   return m;
  176. }
  177.  
  178. gsl_matrix *
  179. create_2x2_matrix(double a11, double a12, double a21, double a22)
  180. {
  181.   gsl_matrix * m = gsl_matrix_alloc(2, 2);
  182.   gsl_matrix_set(m, 0, 0, a11);
  183.   gsl_matrix_set(m, 0, 1, a12);
  184.   gsl_matrix_set(m, 1, 0, a21);
  185.   gsl_matrix_set(m, 1, 1, a22);
  186.   return m;
  187. }
  188.  
  189. gsl_matrix * m35;
  190. gsl_matrix * m53;
  191. gsl_matrix * m97;
  192.  
  193. gsl_matrix * s35;
  194. gsl_matrix * s53;
  195.  
  196. gsl_matrix * hilb2;
  197. gsl_matrix * hilb3;
  198. gsl_matrix * hilb4;
  199. gsl_matrix * hilb12;
  200.  
  201. gsl_matrix * row3;
  202. gsl_matrix * row5;
  203. gsl_matrix * row12;
  204.  
  205. gsl_matrix * A22;
  206.  
  207. gsl_matrix_complex * c7;
  208.  
  209. double m53_lssolution[] = {52.5992295702070, -337.7263113752073, 
  210.                            351.8823436427604};
  211. double hilb2_solution[] = {-8.0, 18.0} ;
  212. double hilb3_solution[] = {27.0, -192.0, 210.0};
  213. double hilb4_solution[] = {-64.0, 900.0, -2520.0, 1820.0};
  214. double hilb12_solution[] = {-1728.0, 245388.0, -8528520.0, 
  215.                             127026900.0, -1009008000.0, 4768571808.0, 
  216.                             -14202796608.0, 27336497760.0, -33921201600.0,
  217.                             26189163000.0, -11437874448.0, 2157916488.0 };
  218.  
  219. double c7_solution[] = { 2.40717272023734e+01, -9.84612797621247e+00,
  220.                          -2.69338853034031e+02, 8.75455232472528e+01,
  221.                          2.96661356736296e+03, -1.02624473923993e+03,
  222.                          -1.82073812124749e+04, 5.67384473042410e+03,
  223.                          5.57693879019068e+04, -1.61540963210502e+04,
  224.                          -7.88941207561151e+04, 1.95053812987858e+04,
  225.                          3.95548551241728e+04, -7.76593696255317e+03 };
  226.  
  227. gsl_matrix * vander2;
  228. gsl_matrix * vander3;
  229. gsl_matrix * vander4;
  230. gsl_matrix * vander12;
  231.  
  232. double vander2_solution[] = {1.0, 0.0}; 
  233. double vander3_solution[] = {0.0, 1.0, 0.0}; 
  234. double vander4_solution[] = {0.0, 0.0, 1.0, 0.0}; 
  235. double vander12_solution[] = {0.0, 0.0, 0.0, 0.0,
  236.                             0.0, 0.0, 0.0, 0.0, 
  237.                             0.0, 0.0, 1.0, 0.0}; 
  238.  
  239. gsl_matrix * moler10;
  240.  
  241. /* matmult now obsolete */
  242. #ifdef MATMULT
  243. int
  244. test_matmult(void)
  245. {
  246.   int s = 0;
  247.  
  248.   gsl_matrix * A = gsl_matrix_calloc(2, 2);
  249.   gsl_matrix * B = gsl_matrix_calloc(2, 3);
  250.   gsl_matrix * C = gsl_matrix_calloc(2, 3);
  251.  
  252.   gsl_matrix_set(A, 0, 0, 10.0);
  253.   gsl_matrix_set(A, 0, 1,  5.0);
  254.   gsl_matrix_set(A, 1, 0,  1.0);
  255.   gsl_matrix_set(A, 1, 1, 20.0);
  256.  
  257.   gsl_matrix_set(B, 0, 0, 10.0);
  258.   gsl_matrix_set(B, 0, 1,  5.0);
  259.   gsl_matrix_set(B, 0, 2,  2.0);
  260.   gsl_matrix_set(B, 1, 0,  1.0);
  261.   gsl_matrix_set(B, 1, 1,  3.0);
  262.   gsl_matrix_set(B, 1, 2,  2.0);
  263.  
  264.   gsl_linalg_matmult(A, B, C);
  265.  
  266.   s += ( fabs(gsl_matrix_get(C, 0, 0) - 105.0) > GSL_DBL_EPSILON );
  267.   s += ( fabs(gsl_matrix_get(C, 0, 1) -  65.0) > GSL_DBL_EPSILON );
  268.   s += ( fabs(gsl_matrix_get(C, 0, 2) -  30.0) > GSL_DBL_EPSILON );
  269.   s += ( fabs(gsl_matrix_get(C, 1, 0) -  30.0) > GSL_DBL_EPSILON );
  270.   s += ( fabs(gsl_matrix_get(C, 1, 1) -  65.0) > GSL_DBL_EPSILON );
  271.   s += ( fabs(gsl_matrix_get(C, 1, 2) -  42.0) > GSL_DBL_EPSILON );
  272.  
  273.   gsl_matrix_free(A);
  274.   gsl_matrix_free(B);
  275.   gsl_matrix_free(C);
  276.  
  277.   return s;
  278. }
  279.  
  280.  
  281. int
  282. test_matmult_mod(void)
  283. {
  284.   int s = 0;
  285.  
  286.   gsl_matrix * A = gsl_matrix_calloc(3, 3);
  287.   gsl_matrix * B = gsl_matrix_calloc(3, 3);
  288.   gsl_matrix * C = gsl_matrix_calloc(3, 3);
  289.   gsl_matrix * D = gsl_matrix_calloc(2, 3);
  290.   gsl_matrix * E = gsl_matrix_calloc(2, 3);
  291.   gsl_matrix * F = gsl_matrix_calloc(2, 2);
  292.  
  293.   gsl_matrix_set(A, 0, 0, 10.0);
  294.   gsl_matrix_set(A, 0, 1,  5.0);
  295.   gsl_matrix_set(A, 0, 2,  1.0);
  296.   gsl_matrix_set(A, 1, 0,  1.0);
  297.   gsl_matrix_set(A, 1, 1, 20.0);
  298.   gsl_matrix_set(A, 1, 2,  5.0);
  299.   gsl_matrix_set(A, 2, 0,  1.0);
  300.   gsl_matrix_set(A, 2, 1,  3.0);
  301.   gsl_matrix_set(A, 2, 2,  7.0);
  302.  
  303.   gsl_matrix_set(B, 0, 0, 10.0);
  304.   gsl_matrix_set(B, 0, 1,  5.0);
  305.   gsl_matrix_set(B, 0, 2,  2.0);
  306.   gsl_matrix_set(B, 1, 0,  1.0);
  307.   gsl_matrix_set(B, 1, 1,  3.0);
  308.   gsl_matrix_set(B, 1, 2,  2.0);
  309.   gsl_matrix_set(B, 2, 0,  1.0);
  310.   gsl_matrix_set(B, 2, 1,  3.0);
  311.   gsl_matrix_set(B, 2, 2,  2.0);
  312.  
  313.   gsl_matrix_set(D, 0, 0, 10.0);
  314.   gsl_matrix_set(D, 0, 1,  5.0);
  315.   gsl_matrix_set(D, 0, 2,  1.0);
  316.   gsl_matrix_set(D, 1, 0,  1.0);
  317.   gsl_matrix_set(D, 1, 1, 20.0);
  318.   gsl_matrix_set(D, 1, 2,  5.0);
  319.  
  320.   gsl_matrix_set(E, 0, 0, 10.0);
  321.   gsl_matrix_set(E, 0, 1,  5.0);
  322.   gsl_matrix_set(E, 0, 2,  2.0);
  323.   gsl_matrix_set(E, 1, 0,  1.0);
  324.   gsl_matrix_set(E, 1, 1,  3.0);
  325.   gsl_matrix_set(E, 1, 2,  2.0);
  326.  
  327.   gsl_linalg_matmult_mod(A, GSL_LINALG_MOD_NONE, B, GSL_LINALG_MOD_NONE, C);
  328.   s += ( fabs(gsl_matrix_get(C, 0, 0) - 106.0) > GSL_DBL_EPSILON );
  329.   s += ( fabs(gsl_matrix_get(C, 0, 1) -  68.0) > GSL_DBL_EPSILON );
  330.   s += ( fabs(gsl_matrix_get(C, 0, 2) -  32.0) > GSL_DBL_EPSILON );
  331.   s += ( fabs(gsl_matrix_get(C, 1, 0) -  35.0) > GSL_DBL_EPSILON );
  332.   s += ( fabs(gsl_matrix_get(C, 1, 1) -  80.0) > GSL_DBL_EPSILON );
  333.   s += ( fabs(gsl_matrix_get(C, 1, 2) -  52.0) > GSL_DBL_EPSILON );
  334.   s += ( fabs(gsl_matrix_get(C, 2, 0) -  20.0) > GSL_DBL_EPSILON );
  335.   s += ( fabs(gsl_matrix_get(C, 2, 1) -  35.0) > GSL_DBL_EPSILON );
  336.   s += ( fabs(gsl_matrix_get(C, 2, 2) -  22.0) > GSL_DBL_EPSILON );
  337.  
  338.   gsl_linalg_matmult_mod(A, GSL_LINALG_MOD_TRANSPOSE, B, GSL_LINALG_MOD_NONE, C);
  339.   s += ( fabs(gsl_matrix_get(C, 0, 0) - 102.0) > GSL_DBL_EPSILON );
  340.   s += ( fabs(gsl_matrix_get(C, 0, 1) -  56.0) > GSL_DBL_EPSILON );
  341.   s += ( fabs(gsl_matrix_get(C, 0, 2) -  24.0) > GSL_DBL_EPSILON );
  342.   s += ( fabs(gsl_matrix_get(C, 1, 0) -  73.0) > GSL_DBL_EPSILON );
  343.   s += ( fabs(gsl_matrix_get(C, 1, 1) -  94.0) > GSL_DBL_EPSILON );
  344.   s += ( fabs(gsl_matrix_get(C, 1, 2) -  56.0) > GSL_DBL_EPSILON );
  345.   s += ( fabs(gsl_matrix_get(C, 2, 0) -  22.0) > GSL_DBL_EPSILON );
  346.   s += ( fabs(gsl_matrix_get(C, 2, 1) -  41.0) > GSL_DBL_EPSILON );
  347.   s += ( fabs(gsl_matrix_get(C, 2, 2) -  26.0) > GSL_DBL_EPSILON );
  348.  
  349.   gsl_linalg_matmult_mod(A, GSL_LINALG_MOD_NONE, B, GSL_LINALG_MOD_TRANSPOSE, C);
  350.   s += ( fabs(gsl_matrix_get(C, 0, 0) - 127.0) > GSL_DBL_EPSILON );
  351.   s += ( fabs(gsl_matrix_get(C, 0, 1) -  27.0) > GSL_DBL_EPSILON );
  352.   s += ( fabs(gsl_matrix_get(C, 0, 2) -  27.0) > GSL_DBL_EPSILON );
  353.   s += ( fabs(gsl_matrix_get(C, 1, 0) - 120.0) > GSL_DBL_EPSILON );
  354.   s += ( fabs(gsl_matrix_get(C, 1, 1) -  71.0) > GSL_DBL_EPSILON );
  355.   s += ( fabs(gsl_matrix_get(C, 1, 2) -  71.0) > GSL_DBL_EPSILON );
  356.   s += ( fabs(gsl_matrix_get(C, 2, 0) -  39.0) > GSL_DBL_EPSILON );
  357.   s += ( fabs(gsl_matrix_get(C, 2, 1) -  24.0) > GSL_DBL_EPSILON );
  358.   s += ( fabs(gsl_matrix_get(C, 2, 2) -  24.0) > GSL_DBL_EPSILON );
  359.  
  360.   gsl_linalg_matmult_mod(A, GSL_LINALG_MOD_TRANSPOSE, B, GSL_LINALG_MOD_TRANSPOSE, C);
  361.   s += ( fabs(gsl_matrix_get(C, 0, 0) - 107.0) > GSL_DBL_EPSILON );
  362.   s += ( fabs(gsl_matrix_get(C, 0, 1) -  15.0) > GSL_DBL_EPSILON );
  363.   s += ( fabs(gsl_matrix_get(C, 0, 2) -  15.0) > GSL_DBL_EPSILON );
  364.   s += ( fabs(gsl_matrix_get(C, 1, 0) - 156.0) > GSL_DBL_EPSILON );
  365.   s += ( fabs(gsl_matrix_get(C, 1, 1) -  71.0) > GSL_DBL_EPSILON );
  366.   s += ( fabs(gsl_matrix_get(C, 1, 2) -  71.0) > GSL_DBL_EPSILON );
  367.   s += ( fabs(gsl_matrix_get(C, 2, 0) -  49.0) > GSL_DBL_EPSILON );
  368.   s += ( fabs(gsl_matrix_get(C, 2, 1) -  30.0) > GSL_DBL_EPSILON );
  369.   s += ( fabs(gsl_matrix_get(C, 2, 2) -  30.0) > GSL_DBL_EPSILON );
  370.  
  371.   /* now try for non-symmetric matrices */
  372.   gsl_linalg_matmult_mod(D, GSL_LINALG_MOD_TRANSPOSE, E, GSL_LINALG_MOD_NONE, C);
  373.   s += ( fabs(gsl_matrix_get(C, 0, 0) - 101.0) > GSL_DBL_EPSILON );
  374.   s += ( fabs(gsl_matrix_get(C, 0, 1) -  53.0) > GSL_DBL_EPSILON );
  375.   s += ( fabs(gsl_matrix_get(C, 0, 2) -  22.0) > GSL_DBL_EPSILON );
  376.   s += ( fabs(gsl_matrix_get(C, 1, 0) -  70.0) > GSL_DBL_EPSILON );
  377.   s += ( fabs(gsl_matrix_get(C, 1, 1) -  85.0) > GSL_DBL_EPSILON );
  378.   s += ( fabs(gsl_matrix_get(C, 1, 2) -  50.0) > GSL_DBL_EPSILON );
  379.   s += ( fabs(gsl_matrix_get(C, 2, 0) -  15.0) > GSL_DBL_EPSILON );
  380.   s += ( fabs(gsl_matrix_get(C, 2, 1) -  20.0) > GSL_DBL_EPSILON );
  381.   s += ( fabs(gsl_matrix_get(C, 2, 2) -  12.0) > GSL_DBL_EPSILON );
  382.  
  383.  
  384.   gsl_linalg_matmult_mod(D, GSL_LINALG_MOD_NONE, E, GSL_LINALG_MOD_TRANSPOSE, F);
  385.   s += ( fabs(gsl_matrix_get(F, 0, 0) - 127.0) > GSL_DBL_EPSILON );
  386.   s += ( fabs(gsl_matrix_get(F, 0, 1) -  27.0) > GSL_DBL_EPSILON );
  387.   s += ( fabs(gsl_matrix_get(F, 1, 0) - 120.0) > GSL_DBL_EPSILON );
  388.   s += ( fabs(gsl_matrix_get(F, 1, 1) -  71.0) > GSL_DBL_EPSILON );
  389.  
  390.  
  391.   gsl_matrix_free(A);
  392.   gsl_matrix_free(B);
  393.   gsl_matrix_free(C);
  394.   gsl_matrix_free(D);
  395.   gsl_matrix_free(E);
  396.   gsl_matrix_free(F);
  397.  
  398.   return s;
  399. }
  400. #endif
  401.  
  402. int
  403. test_LU_solve_dim(const gsl_matrix * m, const double * actual, double eps)
  404. {
  405.   int s = 0;
  406.   int signum;
  407.   size_t i, dim = m->size1;
  408.  
  409.   gsl_permutation * perm = gsl_permutation_alloc(dim);
  410.   gsl_vector * rhs = gsl_vector_alloc(dim);
  411.   gsl_matrix * lu  = gsl_matrix_alloc(dim,dim);
  412.   gsl_vector * x = gsl_vector_alloc(dim);
  413.   gsl_vector * residual = gsl_vector_alloc(dim);
  414.   gsl_matrix_memcpy(lu,m);
  415.   for(i=0; i<dim; i++) gsl_vector_set(rhs, i, i+1.0);
  416.   s += gsl_linalg_LU_decomp(lu, perm, &signum);
  417.   s += gsl_linalg_LU_solve(lu, perm, rhs, x);
  418.  
  419.   for(i=0; i<dim; i++) {
  420.     int foo = check(gsl_vector_get(x, i),actual[i],eps);
  421.     if(foo) {
  422.       printf("%3d[%d]: %22.18g   %22.18g\n", dim, i, gsl_vector_get(x, i), actual[i]);
  423.     }
  424.     s += foo;
  425.   }
  426.  
  427.   s += gsl_linalg_LU_refine(m, lu, perm, rhs, x, residual);
  428.  
  429.   for(i=0; i<dim; i++) {
  430.     int foo = check(gsl_vector_get(x, i),actual[i],eps);
  431.     if(foo) {
  432.       printf("%3d[%d]: %22.18g   %22.18g (improved)\n", dim, i, gsl_vector_get(x, i), actual[i]);
  433.     }
  434.     s += foo;
  435.   }
  436.  
  437.   gsl_vector_free(residual);
  438.   gsl_vector_free(x);
  439.   gsl_matrix_free(lu);
  440.   gsl_vector_free(rhs);
  441.   gsl_permutation_free(perm);
  442.  
  443.   return s;
  444. }
  445.  
  446.  
  447. int test_LU_solve(void)
  448. {
  449.   int f;
  450.   int s = 0;
  451.  
  452.   f = test_LU_solve_dim(hilb2, hilb2_solution, 8.0 * GSL_DBL_EPSILON);
  453.   gsl_test(f, "  LU_solve hilbert(2)");
  454.   s += f;
  455.  
  456.   f = test_LU_solve_dim(hilb3, hilb3_solution, 64.0 * GSL_DBL_EPSILON);
  457.   gsl_test(f, "  LU_solve hilbert(3)");
  458.   s += f;
  459.  
  460.   f = test_LU_solve_dim(hilb4, hilb4_solution, 1024.0 * GSL_DBL_EPSILON);
  461.   gsl_test(f, "  LU_solve hilbert(4)");
  462.   s += f;
  463.  
  464.   f = test_LU_solve_dim(hilb12, hilb12_solution, 0.5);
  465.   gsl_test(f, "  LU_solve hilbert(12)");
  466.   s += f;
  467.  
  468.   f = test_LU_solve_dim(vander2, vander2_solution, 8.0 * GSL_DBL_EPSILON);
  469.   gsl_test(f, "  LU_solve vander(2)");
  470.   s += f;
  471.  
  472.   f = test_LU_solve_dim(vander3, vander3_solution, 64.0 * GSL_DBL_EPSILON);
  473.   gsl_test(f, "  LU_solve vander(3)");
  474.   s += f;
  475.  
  476.   f = test_LU_solve_dim(vander4, vander4_solution, 1024.0 * GSL_DBL_EPSILON);
  477.   gsl_test(f, "  LU_solve vander(4)");
  478.   s += f;
  479.  
  480.   f = test_LU_solve_dim(vander12, vander12_solution, 0.05);
  481.   gsl_test(f, "  LU_solve vander(12)");
  482.   s += f;
  483.  
  484.   return s;
  485. }
  486.  
  487.  
  488. int
  489. test_LUc_solve_dim(const gsl_matrix_complex * m, const double * actual, double eps)
  490. {
  491.   int s = 0;
  492.   int signum;
  493.   size_t i, dim = m->size1;
  494.  
  495.   gsl_permutation * perm = gsl_permutation_alloc(dim);
  496.   gsl_vector_complex * rhs = gsl_vector_complex_alloc(dim);
  497.   gsl_matrix_complex * lu  = gsl_matrix_complex_alloc(dim,dim);
  498.   gsl_vector_complex * x = gsl_vector_complex_alloc(dim);
  499.   gsl_vector_complex * residual = gsl_vector_complex_alloc(dim);
  500.   gsl_matrix_complex_memcpy(lu,m);
  501.   for(i=0; i<dim; i++) 
  502.     {
  503.       gsl_complex z = gsl_complex_rect (2.0*i+1.0, 2.0*i+2.0);
  504.       gsl_vector_complex_set(rhs, i, z);
  505.     }
  506.   s += gsl_linalg_complex_LU_decomp(lu, perm, &signum);
  507.   s += gsl_linalg_complex_LU_solve(lu, perm, rhs, x);
  508.  
  509.   for(i=0; i<dim; i++) {
  510.     gsl_complex z = gsl_vector_complex_get(x, i);
  511.     int foo_r = check(GSL_REAL(z),actual[2*i],eps);
  512.     int foo_i = check(GSL_IMAG(z),actual[2*i+1],eps);
  513.     if(foo_r || foo_i) {
  514.       printf("%3d[%d]: %22.18g   %22.18g\n", dim, i, GSL_REAL(z), actual[2*i]);
  515.       printf("%3d[%d]: %22.18g   %22.18g\n", dim, i, GSL_IMAG(z), actual[2*i+1]);
  516.     }
  517.     s += foo_r + foo_i;
  518.   }
  519.  
  520.   s += gsl_linalg_complex_LU_refine(m, lu, perm, rhs, x, residual);
  521.  
  522.   for(i=0; i<dim; i++) {
  523.     gsl_complex z = gsl_vector_complex_get(x, i);
  524.     int foo_r = check(GSL_REAL(z),actual[2*i],eps);
  525.     int foo_i = check(GSL_IMAG(z),actual[2*i+1],eps);
  526.     if(foo_r || foo_i) {
  527.       printf("%3d[%d]: %22.18g   %22.18g (improved)\n", dim, i, GSL_REAL(z), actual[2*i]);
  528.       printf("%3d[%d]: %22.18g   %22.18g (improved)\n", dim, i, GSL_IMAG(z), actual[2*i+1]);
  529.     }
  530.     s += foo_r + foo_i;
  531.   }
  532.  
  533.   gsl_vector_complex_free(residual);
  534.   gsl_vector_complex_free(x);
  535.   gsl_matrix_complex_free(lu);
  536.   gsl_vector_complex_free(rhs);
  537.   gsl_permutation_free(perm);
  538.  
  539.   return s;
  540. }
  541.  
  542.  
  543. int test_LUc_solve(void)
  544. {
  545.   int f;
  546.   int s = 0;
  547.  
  548.   f = test_LUc_solve_dim(c7, c7_solution, 1024.0 * 1024.0 * GSL_DBL_EPSILON);
  549.   gsl_test(f, "  complex_LU_solve complex(7)");
  550.   s += f;
  551.  
  552.   return s;
  553. }
  554.  
  555.  
  556. int
  557. test_QR_solve_dim(const gsl_matrix * m, const double * actual, double eps)
  558. {
  559.   int s = 0;
  560.   size_t i, dim = m->size1;
  561.  
  562.   gsl_vector * rhs = gsl_vector_alloc(dim);
  563.   gsl_matrix * qr  = gsl_matrix_alloc(dim,dim);
  564.   gsl_vector * d = gsl_vector_alloc(dim);
  565.   gsl_vector * x = gsl_vector_alloc(dim);
  566.  
  567.   gsl_matrix_memcpy(qr,m);
  568.   for(i=0; i<dim; i++) gsl_vector_set(rhs, i, i+1.0);
  569.   s += gsl_linalg_QR_decomp(qr, d);
  570.   s += gsl_linalg_QR_solve(qr, d, rhs, x);
  571.   for(i=0; i<dim; i++) {
  572.     int foo = check(gsl_vector_get(x, i), actual[i], eps);
  573.     if(foo) {
  574.       printf("%3d[%d]: %22.18g   %22.18g\n", dim, i, gsl_vector_get(x, i), actual[i]);
  575.     }
  576.     s += foo;
  577.   }
  578.  
  579.   gsl_vector_free(x);
  580.   gsl_vector_free(d);
  581.   gsl_matrix_free(qr);
  582.   gsl_vector_free(rhs);
  583.  
  584.   return s;
  585. }
  586.  
  587. int test_QR_solve(void)
  588. {
  589.   int f;
  590.   int s = 0;
  591.  
  592.   f = test_QR_solve_dim(hilb2, hilb2_solution, 2 * 8.0 * GSL_DBL_EPSILON);
  593.   gsl_test(f, "  QR_solve hilbert(2)");
  594.   s += f;
  595.  
  596.   f = test_QR_solve_dim(hilb3, hilb3_solution, 2 * 64.0 * GSL_DBL_EPSILON);
  597.   gsl_test(f, "  QR_solve hilbert(3)");
  598.   s += f;
  599.  
  600.   f = test_QR_solve_dim(hilb4, hilb4_solution, 2 * 1024.0 * GSL_DBL_EPSILON);
  601.   gsl_test(f, "  QR_solve hilbert(4)");
  602.   s += f;
  603.  
  604.   f = test_QR_solve_dim(hilb12, hilb12_solution, 0.5);
  605.   gsl_test(f, "  QR_solve hilbert(12)");
  606.   s += f;
  607.  
  608.   f = test_QR_solve_dim(vander2, vander2_solution, 8.0 * GSL_DBL_EPSILON);
  609.   gsl_test(f, "  QR_solve vander(2)");
  610.   s += f;
  611.  
  612.   f = test_QR_solve_dim(vander3, vander3_solution, 64.0 * GSL_DBL_EPSILON);
  613.   gsl_test(f, "  QR_solve vander(3)");
  614.   s += f;
  615.  
  616.   f = test_QR_solve_dim(vander4, vander4_solution, 1024.0 * GSL_DBL_EPSILON);
  617.   gsl_test(f, "  QR_solve vander(4)");
  618.   s += f;
  619.  
  620.   f = test_QR_solve_dim(vander12, vander12_solution, 0.05);
  621.   gsl_test(f, "  QR_solve vander(12)");
  622.   s += f;
  623.  
  624.   return s;
  625. }
  626.  
  627.  
  628. int
  629. test_QR_lssolve_dim(const gsl_matrix * m, const double * actual, double eps)
  630. {
  631.   int s = 0;
  632.   size_t i, M = m->size1, N = m->size2;
  633.  
  634.   gsl_vector * rhs = gsl_vector_alloc(M);
  635.   gsl_matrix * qr  = gsl_matrix_alloc(M,N);
  636.   gsl_vector * d = gsl_vector_alloc(N);
  637.   gsl_vector * x = gsl_vector_alloc(N);
  638.   gsl_vector * r = gsl_vector_alloc(M);
  639.   gsl_vector * res = gsl_vector_alloc(M);
  640.  
  641.   gsl_matrix_memcpy(qr,m);
  642.   for(i=0; i<M; i++) gsl_vector_set(rhs, i, i+1.0);
  643.   s += gsl_linalg_QR_decomp(qr, d);
  644.   s += gsl_linalg_QR_lssolve(qr, d, rhs, x, res);
  645.  
  646.   for(i=0; i<N; i++) {
  647.     int foo = check(gsl_vector_get(x, i), actual[i], eps);
  648.     if(foo) {
  649.       printf("(%3d,%3d)[%d]: %22.18g   %22.18g\n", M, N, i, gsl_vector_get(x, i), actual[i]);
  650.     }
  651.     s += foo;
  652.   }
  653.  
  654.   /* compute residual r = b - m x */
  655.   if (M == N) {
  656.     gsl_vector_set_zero(r);
  657.   } else {
  658.     gsl_vector_memcpy(r, rhs);
  659.     gsl_blas_dgemv(CblasNoTrans, -1.0, m, x, 1.0, r);
  660.   };
  661.  
  662.   for(i=0; i<N; i++) {
  663.     int foo = check(gsl_vector_get(res, i), gsl_vector_get(r,i), sqrt(eps));
  664.     if(foo) {
  665.       printf("(%3d,%3d)[%d]: %22.18g   %22.18g\n", M, N, i, gsl_vector_get(res, i), gsl_vector_get(r,i));
  666.     }
  667.     s += foo;
  668.   }
  669.  
  670.   gsl_vector_free(r);
  671.   gsl_vector_free(res);
  672.   gsl_vector_free(x);
  673.   gsl_vector_free(d);
  674.   gsl_matrix_free(qr);
  675.   gsl_vector_free(rhs);
  676.  
  677.   return s;
  678. }
  679.  
  680. int test_QR_lssolve(void)
  681. {
  682.   int f;
  683.   int s = 0;
  684.  
  685.   f = test_QR_lssolve_dim(m53, m53_lssolution, 2 * 64.0 * GSL_DBL_EPSILON);
  686.   gsl_test(f, "  QR_lssolve m(5,3)");
  687.   s += f;
  688.  
  689.   f = test_QR_lssolve_dim(hilb2, hilb2_solution, 2 * 8.0 * GSL_DBL_EPSILON);
  690.   gsl_test(f, "  QR_lssolve hilbert(2)");
  691.   s += f;
  692.  
  693.   f = test_QR_lssolve_dim(hilb3, hilb3_solution, 2 * 64.0 * GSL_DBL_EPSILON);
  694.   gsl_test(f, "  QR_lssolve hilbert(3)");
  695.   s += f;
  696.  
  697.   f = test_QR_lssolve_dim(hilb4, hilb4_solution, 2 * 1024.0 * GSL_DBL_EPSILON);
  698.   gsl_test(f, "  QR_lssolve hilbert(4)");
  699.   s += f;
  700.  
  701.   f = test_QR_lssolve_dim(hilb12, hilb12_solution, 0.5);
  702.   gsl_test(f, "  QR_lssolve hilbert(12)");
  703.   s += f;
  704.  
  705.   f = test_QR_lssolve_dim(vander2, vander2_solution, 8.0 * GSL_DBL_EPSILON);
  706.   gsl_test(f, "  QR_lssolve vander(2)");
  707.   s += f;
  708.  
  709.   f = test_QR_lssolve_dim(vander3, vander3_solution, 64.0 * GSL_DBL_EPSILON);
  710.   gsl_test(f, "  QR_lssolve vander(3)");
  711.   s += f;
  712.  
  713.   f = test_QR_lssolve_dim(vander4, vander4_solution, 1024.0 * GSL_DBL_EPSILON);
  714.   gsl_test(f, "  QR_lssolve vander(4)");
  715.   s += f;
  716.  
  717.   f = test_QR_lssolve_dim(vander12, vander12_solution, 0.05);
  718.   gsl_test(f, "  QR_lssolve vander(12)");
  719.   s += f;
  720.  
  721.   return s;
  722. }
  723.  
  724.  
  725. int
  726. test_QR_decomp_dim(const gsl_matrix * m, double eps)
  727. {
  728.   int s = 0;
  729.   size_t i,j, M = m->size1, N = m->size2;
  730.  
  731.   gsl_matrix * qr = gsl_matrix_alloc(M,N);
  732.   gsl_matrix * a  = gsl_matrix_alloc(M,N);
  733.   gsl_matrix * q  = gsl_matrix_alloc(M,M);
  734.   gsl_matrix * r  = gsl_matrix_alloc(M,N);
  735.   gsl_vector * d = gsl_vector_alloc(GSL_MIN(M,N));
  736.  
  737.   gsl_matrix_memcpy(qr,m);
  738.  
  739.   s += gsl_linalg_QR_decomp(qr, d);
  740.   s += gsl_linalg_QR_unpack(qr, d, q, r);
  741.   
  742.   /* compute a = q r */
  743.   gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, q, r, 0.0, a);
  744.  
  745.   for(i=0; i<M; i++) {
  746.     for(j=0; j<N; j++) {
  747.       double aij = gsl_matrix_get(a, i, j);
  748.       double mij = gsl_matrix_get(m, i, j);
  749.       int foo = check(aij, mij, eps);
  750.       if(foo) {
  751.         printf("(%3d,%3d)[%d,%d]: %22.18g   %22.18g\n", M, N, i,j, aij, mij);
  752.       }
  753.       s += foo;
  754.     }
  755.   }
  756.  
  757.   gsl_vector_free(d);
  758.   gsl_matrix_free(qr);
  759.   gsl_matrix_free(a);
  760.   gsl_matrix_free(q);
  761.   gsl_matrix_free(r);
  762.  
  763.   return s;
  764. }
  765.  
  766. int test_QR_decomp(void)
  767. {
  768.   int f;
  769.   int s = 0;
  770.  
  771.   f = test_QR_decomp_dim(m35, 2 * 8.0 * GSL_DBL_EPSILON);
  772.   gsl_test(f, "  QR_decomp m(3,5)");
  773.   s += f;
  774.  
  775.   f = test_QR_decomp_dim(m53, 2 * 64.0 * GSL_DBL_EPSILON);
  776.   gsl_test(f, "  QR_decomp m(5,3)");
  777.   s += f;
  778.  
  779.   f = test_QR_decomp_dim(hilb2, 2 * 8.0 * GSL_DBL_EPSILON);
  780.   gsl_test(f, "  QR_decomp hilbert(2)");
  781.   s += f;
  782.  
  783.   f = test_QR_decomp_dim(hilb3, 2 * 64.0 * GSL_DBL_EPSILON);
  784.   gsl_test(f, "  QR_decomp hilbert(3)");
  785.   s += f;
  786.  
  787.   f = test_QR_decomp_dim(hilb4, 2 * 1024.0 * GSL_DBL_EPSILON);
  788.   gsl_test(f, "  QR_decomp hilbert(4)");
  789.   s += f;
  790.  
  791.   f = test_QR_decomp_dim(hilb12, 2 * 1024.0 * GSL_DBL_EPSILON);
  792.   gsl_test(f, "  QR_decomp hilbert(12)");
  793.   s += f;
  794.  
  795.   f = test_QR_decomp_dim(vander2, 8.0 * GSL_DBL_EPSILON);
  796.   gsl_test(f, "  QR_decomp vander(2)");
  797.   s += f;
  798.  
  799.   f = test_QR_decomp_dim(vander3, 64.0 * GSL_DBL_EPSILON);
  800.   gsl_test(f, "  QR_decomp vander(3)");
  801.   s += f;
  802.  
  803.   f = test_QR_decomp_dim(vander4, 1024.0 * GSL_DBL_EPSILON);
  804.   gsl_test(f, "  QR_decomp vander(4)");
  805.   s += f;
  806.  
  807.   f = test_QR_decomp_dim(vander12, 0.0005); /* FIXME: bad accuracy */
  808.   gsl_test(f, "  QR_decomp vander(12)");
  809.   s += f;
  810.  
  811.   return s;
  812. }
  813.  
  814. int
  815. test_QRPT_solve_dim(const gsl_matrix * m, const double * actual, double eps)
  816. {
  817.   int s = 0;
  818.   int signum;
  819.   size_t i, dim = m->size1;
  820.  
  821.   gsl_permutation * perm = gsl_permutation_alloc(dim);
  822.   gsl_vector * rhs = gsl_vector_alloc(dim);
  823.   gsl_matrix * qr  = gsl_matrix_alloc(dim,dim);
  824.   gsl_vector * d = gsl_vector_alloc(dim);
  825.   gsl_vector * x = gsl_vector_alloc(dim);
  826.   gsl_vector * norm = gsl_vector_alloc(dim);
  827.  
  828.   gsl_matrix_memcpy(qr,m);
  829.   for(i=0; i<dim; i++) gsl_vector_set(rhs, i, i+1.0);
  830.   s += gsl_linalg_QRPT_decomp(qr, d, perm, &signum, norm);
  831.   s += gsl_linalg_QRPT_solve(qr, d, perm, rhs, x);
  832.   for(i=0; i<dim; i++) {
  833.     int foo = check(gsl_vector_get(x, i), actual[i], eps);
  834.     if(foo) {
  835.       printf("%3d[%d]: %22.18g   %22.18g\n", dim, i, gsl_vector_get(x, i), actual[i]);
  836.     }
  837.     s += foo;
  838.   }
  839.  
  840.   gsl_vector_free(norm);
  841.   gsl_vector_free(x);
  842.   gsl_vector_free(d);
  843.   gsl_matrix_free(qr);
  844.   gsl_vector_free(rhs);
  845.   gsl_permutation_free(perm);
  846.  
  847.   return s;
  848. }
  849.  
  850. int test_QRPT_solve(void)
  851. {
  852.   int f;
  853.   int s = 0;
  854.  
  855.   f = test_QRPT_solve_dim(hilb2, hilb2_solution, 2 * 8.0 * GSL_DBL_EPSILON);
  856.   gsl_test(f, "  QRPT_solve hilbert(2)");
  857.   s += f;
  858.  
  859.   f = test_QRPT_solve_dim(hilb3, hilb3_solution, 2 * 64.0 * GSL_DBL_EPSILON);
  860.   gsl_test(f, "  QRPT_solve hilbert(3)");
  861.   s += f;
  862.  
  863.   f = test_QRPT_solve_dim(hilb4, hilb4_solution, 2 * 2048.0 * GSL_DBL_EPSILON);
  864.   gsl_test(f, "  QRPT_solve hilbert(4)");
  865.   s += f;
  866.  
  867.   f = test_QRPT_solve_dim(hilb12, hilb12_solution, 0.5);
  868.   gsl_test(f, "  QRPT_solve hilbert(12)");
  869.   s += f;
  870.  
  871.   f = test_QRPT_solve_dim(vander2, vander2_solution, 8.0 * GSL_DBL_EPSILON);
  872.   gsl_test(f, "  QRPT_solve vander(2)");
  873.   s += f;
  874.  
  875.   f = test_QRPT_solve_dim(vander3, vander3_solution, 64.0 * GSL_DBL_EPSILON);
  876.   gsl_test(f, "  QRPT_solve vander(3)");
  877.   s += f;
  878.  
  879.   f = test_QRPT_solve_dim(vander4, vander4_solution, 1024.0 * GSL_DBL_EPSILON);
  880.   gsl_test(f, "  QRPT_solve vander(4)");
  881.   s += f;
  882.  
  883.   f = test_QRPT_solve_dim(vander12, vander12_solution, 0.05);
  884.   gsl_test(f, "  QRPT_solve vander(12)");
  885.   s += f;
  886.  
  887.   return s;
  888. }
  889.  
  890.  
  891. int
  892. test_QRPT_decomp_dim(const gsl_matrix * m, double eps)
  893. {
  894.   int s = 0, signum;
  895.   size_t i,j, M = m->size1, N = m->size2;
  896.  
  897.   gsl_matrix * qr = gsl_matrix_alloc(M,N);
  898.   gsl_matrix * a  = gsl_matrix_alloc(M,N);
  899.   gsl_matrix * q  = gsl_matrix_alloc(M,M);
  900.   gsl_matrix * r  = gsl_matrix_alloc(M,N);
  901.   gsl_vector * d = gsl_vector_alloc(GSL_MIN(M,N));
  902.   gsl_vector * norm = gsl_vector_alloc(N);
  903.  
  904.   gsl_permutation * perm = gsl_permutation_alloc(N);
  905.  
  906.   gsl_matrix_memcpy(qr,m);
  907.  
  908.   s += gsl_linalg_QRPT_decomp(qr, d, perm, &signum, norm);
  909.   s += gsl_linalg_QR_unpack(qr, d, q, r);
  910.  
  911.   /* compute a = q r */
  912.   gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, q, r, 0.0, a);
  913.  
  914.  
  915.   /* Compute QR P^T by permuting the elements of the rows of QR */
  916.  
  917.   for (i = 0; i < M; i++) {
  918.     gsl_vector_view row = gsl_matrix_row (a, i);
  919.     gsl_permute_vector_inverse (perm, &row.vector);
  920.   }
  921.  
  922.   for(i=0; i<M; i++) {
  923.     for(j=0; j<N; j++) {
  924.       double aij = gsl_matrix_get(a, i, j);
  925.       double mij = gsl_matrix_get(m, i, j);
  926.       int foo = check(aij, mij, eps);
  927.       if(foo) {
  928.         printf("(%3d,%3d)[%d,%d]: %22.18g   %22.18g\n", M, N, i,j, aij, mij);
  929.       }
  930.       s += foo;
  931.     }
  932.   }
  933.  
  934.   gsl_permutation_free (perm);
  935.   gsl_vector_free(norm);
  936.   gsl_vector_free(d);
  937.   gsl_matrix_free(qr);
  938.   gsl_matrix_free(a);
  939.   gsl_matrix_free(q);
  940.   gsl_matrix_free(r);
  941.  
  942.   return s;
  943. }
  944.  
  945. int test_QRPT_decomp(void)
  946. {
  947.   int f;
  948.   int s = 0;
  949.  
  950.   f = test_QRPT_decomp_dim(m35, 2 * 8.0 * GSL_DBL_EPSILON);
  951.   gsl_test(f, "  QRPT_decomp m(3,5)");
  952.   s += f;
  953.  
  954.   f = test_QRPT_decomp_dim(m53, 2 * 8.0 * GSL_DBL_EPSILON);
  955.   gsl_test(f, "  QRPT_decomp m(5,3)");
  956.   s += f;
  957.  
  958.   f = test_QRPT_decomp_dim(s35, 2 * 8.0 * GSL_DBL_EPSILON);
  959.   gsl_test(f, "  QRPT_decomp s(3,5)");
  960.   s += f;
  961.  
  962.   f = test_QRPT_decomp_dim(s53, 2 * 8.0 * GSL_DBL_EPSILON);
  963.   gsl_test(f, "  QRPT_decomp s(5,3)");
  964.   s += f;
  965.  
  966.   f = test_QRPT_decomp_dim(hilb2, 2 * 8.0 * GSL_DBL_EPSILON);
  967.   gsl_test(f, "  QRPT_decomp hilbert(2)");
  968.   s += f;
  969.  
  970.   f = test_QRPT_decomp_dim(hilb3, 2 * 64.0 * GSL_DBL_EPSILON);
  971.   gsl_test(f, "  QRPT_decomp hilbert(3)");
  972.   s += f;
  973.  
  974.   f = test_QRPT_decomp_dim(hilb4, 2 * 1024.0 * GSL_DBL_EPSILON);
  975.   gsl_test(f, "  QRPT_decomp hilbert(4)");
  976.   s += f;
  977.  
  978.   f = test_QRPT_decomp_dim(hilb12, 2 * 1024.0 * GSL_DBL_EPSILON);
  979.   gsl_test(f, "  QRPT_decomp hilbert(12)");
  980.   s += f;
  981.  
  982.   f = test_QRPT_decomp_dim(vander2, 8.0 * GSL_DBL_EPSILON);
  983.   gsl_test(f, "  QRPT_decomp vander(2)");
  984.   s += f;
  985.  
  986.   f = test_QRPT_decomp_dim(vander3, 64.0 * GSL_DBL_EPSILON);
  987.   gsl_test(f, "  QRPT_decomp vander(3)");
  988.   s += f;
  989.  
  990.   f = test_QRPT_decomp_dim(vander4, 1024.0 * GSL_DBL_EPSILON);
  991.   gsl_test(f, "  QRPT_decomp vander(4)");
  992.   s += f;
  993.  
  994.   f = test_QRPT_decomp_dim(vander12, 0.0005); /* FIXME: bad accuracy */
  995.   gsl_test(f, "  QRPT_decomp vander(12)");
  996.   s += f;
  997.  
  998.   return s;
  999. }
  1000.  
  1001.  
  1002. int
  1003. test_QR_update_dim(const gsl_matrix * m, double eps)
  1004. {
  1005.   int s = 0;
  1006.   size_t i,j,k, M = m->size1, N = m->size2;
  1007.  
  1008.   gsl_vector * rhs = gsl_vector_alloc(N);
  1009.   gsl_matrix * qr1  = gsl_matrix_alloc(M,N);
  1010.   gsl_matrix * qr2  = gsl_matrix_alloc(M,N);
  1011.   gsl_matrix * q1  = gsl_matrix_alloc(M,M);
  1012.   gsl_matrix * r1  = gsl_matrix_alloc(M,N);
  1013.   gsl_matrix * q2  = gsl_matrix_alloc(M,M);
  1014.   gsl_matrix * r2  = gsl_matrix_alloc(M,N);
  1015.   gsl_vector * d = gsl_vector_alloc(GSL_MIN(M,N));
  1016.   gsl_vector * solution1 = gsl_vector_alloc(N);
  1017.   gsl_vector * solution2 = gsl_vector_alloc(N);
  1018.   gsl_vector * u = gsl_vector_alloc(M);
  1019.   gsl_vector * v = gsl_vector_alloc(N);
  1020.   gsl_vector * w = gsl_vector_alloc(M);
  1021.  
  1022.   gsl_matrix_memcpy(qr1,m);
  1023.   gsl_matrix_memcpy(qr2,m);
  1024.   for(i=0; i<N; i++) gsl_vector_set(rhs, i, i+1.0);
  1025.   for(i=0; i<M; i++) gsl_vector_set(u, i, sin(i+1.0));
  1026.   for(i=0; i<N; i++) gsl_vector_set(v, i, cos(i+2.0) + sin(i*i+3.0));
  1027.  
  1028.   for(i=0; i<M; i++) 
  1029.     {
  1030.       double ui = gsl_vector_get(u, i);
  1031.       for(j=0; j<N; j++) 
  1032.         {
  1033.           double vj = gsl_vector_get(v, j);
  1034.           double qij = gsl_matrix_get(qr1, i, j);
  1035.           gsl_matrix_set(qr1, i, j, qij + ui * vj);
  1036.         }
  1037.     }
  1038.  
  1039.   s += gsl_linalg_QR_decomp(qr2, d);
  1040.   s += gsl_linalg_QR_unpack(qr2, d, q2, r2);
  1041.  
  1042.   /* compute w = Q^T u */
  1043.       
  1044.   for (j = 0; j < M; j++)
  1045.     {
  1046.       double sum = 0;
  1047.       for (i = 0; i < M; i++)
  1048.           sum += gsl_matrix_get (q2, i, j) * gsl_vector_get (u, i);
  1049.       gsl_vector_set (w, j, sum);
  1050.     }
  1051.  
  1052.   s += gsl_linalg_QR_update(q2, r2, w, v);
  1053.  
  1054.   /* compute qr2 = q2 * r2 */
  1055.  
  1056.   for (i = 0; i < M; i++)
  1057.     {
  1058.       for (j = 0; j< N; j++)
  1059.         {
  1060.           double sum = 0;
  1061.           for (k = 0; k <= GSL_MIN(j,M-1); k++)
  1062.             {
  1063.               double qik = gsl_matrix_get(q2, i, k);
  1064.               double rkj = gsl_matrix_get(r2, k, j);
  1065.               sum += qik * rkj ;
  1066.             }
  1067.           gsl_matrix_set (qr2, i, j, sum);
  1068.         }
  1069.     }
  1070.  
  1071.   for(i=0; i<M; i++) {
  1072.     for(j=0; j<N; j++) {
  1073.       double s1 = gsl_matrix_get(qr1, i, j);
  1074.       double s2 = gsl_matrix_get(qr2, i, j);
  1075.       
  1076.       int foo = check(s1, s2, eps);
  1077.       if(foo) {
  1078.         printf("(%3d,%3d)[%d,%d]: %22.18g   %22.18g\n", M, N, i,j, s1, s2);
  1079.       }
  1080.       s += foo;
  1081.     }
  1082.   }
  1083.  
  1084.   gsl_vector_free(solution1);
  1085.   gsl_vector_free(solution2);
  1086.   gsl_vector_free(d);
  1087.   gsl_vector_free(u);
  1088.   gsl_vector_free(v);
  1089.   gsl_vector_free(w);
  1090.   gsl_matrix_free(qr1);
  1091.   gsl_matrix_free(qr2);
  1092.   gsl_matrix_free(q1);
  1093.   gsl_matrix_free(r1);
  1094.   gsl_matrix_free(q2);
  1095.   gsl_matrix_free(r2);
  1096.   gsl_vector_free(rhs);
  1097.  
  1098.   return s;
  1099. }
  1100.  
  1101. int test_QR_update(void)
  1102. {
  1103.   int f;
  1104.   int s = 0;
  1105.  
  1106.   f = test_QR_update_dim(m35, 2 * 512.0 * GSL_DBL_EPSILON);
  1107.   gsl_test(f, "  QR_update m(3,5)");
  1108.   s += f;
  1109.  
  1110.   f = test_QR_update_dim(m53, 2 * 512.0 * GSL_DBL_EPSILON);
  1111.   gsl_test(f, "  QR_update m(5,3)");
  1112.   s += f;
  1113.  
  1114.   f = test_QR_update_dim(hilb2,  2 * 512.0 * GSL_DBL_EPSILON);
  1115.   gsl_test(f, "  QR_update hilbert(2)");
  1116.   s += f;
  1117.  
  1118.   f = test_QR_update_dim(hilb3,  2 * 512.0 * GSL_DBL_EPSILON);
  1119.   gsl_test(f, "  QR_update hilbert(3)");
  1120.   s += f;
  1121.  
  1122.   f = test_QR_update_dim(hilb4, 2 * 1024.0 * GSL_DBL_EPSILON);
  1123.   gsl_test(f, "  QR_update hilbert(4)");
  1124.   s += f;
  1125.  
  1126.   f = test_QR_update_dim(hilb12, 2 * 1024.0 * GSL_DBL_EPSILON);
  1127.   gsl_test(f, "  QR_update hilbert(12)");
  1128.   s += f;
  1129.  
  1130.   f = test_QR_update_dim(vander2, 8.0 * GSL_DBL_EPSILON);
  1131.   gsl_test(f, "  QR_update vander(2)");
  1132.   s += f;
  1133.  
  1134.   f = test_QR_update_dim(vander3, 64.0 * GSL_DBL_EPSILON);
  1135.   gsl_test(f, "  QR_update vander(3)");
  1136.   s += f;
  1137.  
  1138.   f = test_QR_update_dim(vander4, 1024.0 * GSL_DBL_EPSILON);
  1139.   gsl_test(f, "  QR_update vander(4)");
  1140.   s += f;
  1141.  
  1142.   f = test_QR_update_dim(vander12, 0.0005); /* FIXME: bad accuracy */
  1143.   gsl_test(f, "  QR_update vander(12)");
  1144.   s += f;
  1145.  
  1146.   return s;
  1147. }
  1148.  
  1149.  
  1150. int
  1151. test_SV_solve_dim(const gsl_matrix * m, const double * actual, double eps)
  1152. {
  1153.   int s = 0;
  1154.   size_t i, dim = m->size1;
  1155.  
  1156.   gsl_vector * rhs = gsl_vector_alloc(dim);
  1157.   gsl_matrix * u  = gsl_matrix_alloc(dim,dim);
  1158.   gsl_matrix * q  = gsl_matrix_alloc(dim,dim);
  1159.   gsl_vector * d = gsl_vector_alloc(dim);
  1160.   gsl_vector * x = gsl_vector_calloc(dim);
  1161.   gsl_matrix_memcpy(u,m);
  1162.   for(i=0; i<dim; i++) gsl_vector_set(rhs, i, i+1.0);
  1163.   s += gsl_linalg_SV_decomp(u, q, d, x);
  1164.   s += gsl_linalg_SV_solve(u, q, d, rhs, x);
  1165.   for(i=0; i<dim; i++) {
  1166.     int foo = check(gsl_vector_get(x, i), actual[i], eps);
  1167.     if(foo) {
  1168.       printf("%3d[%d]: %22.18g   %22.18g\n", dim, i, gsl_vector_get(x, i), actual[i]);
  1169.     }
  1170.     s += foo;
  1171.   }
  1172.   gsl_vector_free(x);
  1173.   gsl_vector_free(d);
  1174.   gsl_matrix_free(u);
  1175.   gsl_matrix_free(q);
  1176.   gsl_vector_free(rhs);
  1177.  
  1178.   return s;
  1179. }
  1180.  
  1181. int test_SV_solve(void)
  1182. {
  1183.   int f;
  1184.   int s = 0;
  1185.  
  1186.   f = test_SV_solve_dim(hilb2, hilb2_solution, 2 * 8.0 * GSL_DBL_EPSILON);
  1187.   gsl_test(f, "  SV_solve hilbert(2)");
  1188.   s += f;
  1189.  
  1190.   f = test_SV_solve_dim(hilb3, hilb3_solution, 2 * 64.0 * GSL_DBL_EPSILON);
  1191.   gsl_test(f, "  SV_solve hilbert(3)");
  1192.   s += f;
  1193.  
  1194.   f = test_SV_solve_dim(hilb4, hilb4_solution, 2 * 1024.0 * GSL_DBL_EPSILON);
  1195.   gsl_test(f, "  SV_solve hilbert(4)");
  1196.   s += f;
  1197.  
  1198.   f = test_SV_solve_dim(hilb12, hilb12_solution, 0.5);
  1199.   gsl_test(f, "  SV_solve hilbert(12)");
  1200.   s += f;
  1201.  
  1202.   f = test_SV_solve_dim(vander2, vander2_solution, 64.0 * GSL_DBL_EPSILON);
  1203.   gsl_test(f, "  SV_solve vander(2)");
  1204.   s += f;
  1205.  
  1206.   f = test_SV_solve_dim(vander3, vander3_solution, 64.0 * GSL_DBL_EPSILON);
  1207.   gsl_test(f, "  SV_solve vander(3)");
  1208.   s += f;
  1209.  
  1210.   f = test_SV_solve_dim(vander4, vander4_solution, 1024.0 * GSL_DBL_EPSILON);
  1211.   gsl_test(f, "  SV_solve vander(4)");
  1212.   s += f;
  1213.  
  1214.   f = test_SV_solve_dim(vander12, vander12_solution, 0.05);
  1215.   gsl_test(f, "  SV_solve vander(12)");
  1216.   s += f;
  1217.  
  1218.   return s;
  1219. }
  1220.  
  1221.  
  1222. int
  1223. test_SV_decomp_dim(const gsl_matrix * m, double eps)
  1224. {
  1225.   int s = 0;
  1226.   double di1;
  1227.   size_t i,j, M = m->size1, N = m->size2;
  1228.  
  1229.   gsl_matrix * v  = gsl_matrix_alloc(M,N);
  1230.   gsl_matrix * a  = gsl_matrix_alloc(M,N);
  1231.   gsl_matrix * q  = gsl_matrix_alloc(N,N);
  1232.   gsl_matrix * dqt  = gsl_matrix_alloc(N,N);
  1233.   gsl_vector * d  = gsl_vector_alloc(N);
  1234.   gsl_vector * w  = gsl_vector_alloc(N);
  1235.  
  1236.   gsl_matrix_memcpy(v,m);
  1237.  
  1238.   s += gsl_linalg_SV_decomp(v, q, d, w);
  1239.  
  1240.   /* Check that singular values are non-negative and in non-decreasing
  1241.      order */
  1242.   
  1243.   di1 = 0.0;
  1244.  
  1245.   for (i = 0; i < N; i++)
  1246.     {
  1247.       double di = gsl_vector_get (d, i);
  1248.  
  1249.       if (di < 0) {
  1250.         s++;
  1251.         printf("singular value %d = %22.18g < 0\n", i, di);
  1252.       }
  1253.  
  1254.       if(i > 0 && di > di1) {
  1255.         s++;
  1256.         printf("singular value %d = %22.18g vs previous %22.18g\n", i, di, di1);
  1257.       }
  1258.  
  1259.       di1 = di;
  1260.     }      
  1261.   
  1262.   /* Scale dqt = D Q^T */
  1263.   
  1264.   for (i = 0; i < N ; i++)
  1265.     {
  1266.       double di = gsl_vector_get (d, i);
  1267.  
  1268.       for (j = 0; j < N; j++)
  1269.         {
  1270.           double qji = gsl_matrix_get(q, j, i);
  1271.           gsl_matrix_set (dqt, i, j, qji * di);
  1272.         }
  1273.     }
  1274.             
  1275.   /* compute a = v dqt */
  1276.   gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, v, dqt, 0.0, a);
  1277.  
  1278.   for(i=0; i<M; i++) {
  1279.     for(j=0; j<N; j++) {
  1280.       double aij = gsl_matrix_get(a, i, j);
  1281.       double mij = gsl_matrix_get(m, i, j);
  1282.       int foo = check(aij, mij, eps);
  1283.       if(foo) {
  1284.         printf("(%3d,%3d)[%d,%d]: %22.18g   %22.18g\n", M, N, i,j, aij, mij);
  1285.       }
  1286.       s += foo;
  1287.     }
  1288.   }
  1289.   gsl_vector_free(w);
  1290.   gsl_vector_free(d);
  1291.   gsl_matrix_free(v);
  1292.   gsl_matrix_free(a);
  1293.   gsl_matrix_free(q);
  1294.   gsl_matrix_free(dqt);
  1295.  
  1296.   return s;
  1297. }
  1298.  
  1299. int test_SV_decomp(void)
  1300. {
  1301.   int f;
  1302.   int s = 0;
  1303.  
  1304.   /* M<N not implemented yet */
  1305. #if 0
  1306.   f = test_SV_decomp_dim(m35, 2 * 8.0 * GSL_DBL_EPSILON);
  1307.   gsl_test(f, "  SV_decomp m(3,5)");
  1308.   s += f;
  1309. #endif
  1310.   f = test_SV_decomp_dim(m53, 2 * 64.0 * GSL_DBL_EPSILON);
  1311.   gsl_test(f, "  SV_decomp m(5,3)");
  1312.   s += f;
  1313.  
  1314.   f = test_SV_decomp_dim(moler10, 2 * 64.0 * GSL_DBL_EPSILON);
  1315.   gsl_test(f, "  SV_decomp moler(10)");
  1316.   s += f;
  1317.  
  1318.   f = test_SV_decomp_dim(hilb2, 2 * 8.0 * GSL_DBL_EPSILON);
  1319.   gsl_test(f, "  SV_decomp hilbert(2)");
  1320.   s += f;
  1321.  
  1322.   f = test_SV_decomp_dim(hilb3, 2 * 64.0 * GSL_DBL_EPSILON);
  1323.   gsl_test(f, "  SV_decomp hilbert(3)");
  1324.   s += f;
  1325.  
  1326.   f = test_SV_decomp_dim(hilb4, 2 * 1024.0 * GSL_DBL_EPSILON);
  1327.   gsl_test(f, "  SV_decomp hilbert(4)");
  1328.   s += f;
  1329.  
  1330.   f = test_SV_decomp_dim(hilb12, 2 * 1024.0 * GSL_DBL_EPSILON);
  1331.   gsl_test(f, "  SV_decomp hilbert(12)");
  1332.   s += f;
  1333.  
  1334.   f = test_SV_decomp_dim(vander2, 8.0 * GSL_DBL_EPSILON);
  1335.   gsl_test(f, "  SV_decomp vander(2)");
  1336.   s += f;
  1337.  
  1338.   f = test_SV_decomp_dim(vander3, 64.0 * GSL_DBL_EPSILON);
  1339.   gsl_test(f, "  SV_decomp vander(3)");
  1340.   s += f;
  1341.  
  1342.   f = test_SV_decomp_dim(vander4, 1024.0 * GSL_DBL_EPSILON);
  1343.   gsl_test(f, "  SV_decomp vander(4)");
  1344.   s += f;
  1345.  
  1346.   f = test_SV_decomp_dim(vander12, 1e-4);
  1347.   gsl_test(f, "  SV_decomp vander(12)");
  1348.   s += f;
  1349.  
  1350.   f = test_SV_decomp_dim(row3, 10 * GSL_DBL_EPSILON);
  1351.   gsl_test(f, "  SV_decomp row3");
  1352.   s += f;
  1353.  
  1354.   f = test_SV_decomp_dim(row5, 128 * GSL_DBL_EPSILON);
  1355.   gsl_test(f, "  SV_decomp row5");
  1356.   s += f;
  1357.  
  1358.   f = test_SV_decomp_dim(row12, 1024 * GSL_DBL_EPSILON);
  1359.   gsl_test(f, "  SV_decomp row12");
  1360.   s += f;
  1361.  
  1362.   {
  1363.     double i1, i2, i3, i4;
  1364.     double lower = -2, upper = 2;
  1365.  
  1366.     for (i1 = lower; i1 <= upper; i1++)
  1367.       {
  1368.         for (i2 = lower; i2 <= upper; i2++)
  1369.           {
  1370.             for (i3 = lower; i3 <= upper; i3++)
  1371.               {
  1372.                 for (i4 = lower; i4 <= upper; i4++)
  1373.                   {
  1374.                     gsl_matrix_set (A22, 0,0, i1);
  1375.                     gsl_matrix_set (A22, 0,1, i2);
  1376.                     gsl_matrix_set (A22, 1,0, i3);
  1377.                     gsl_matrix_set (A22, 1,1, i4);
  1378.                     
  1379.                     f = test_SV_decomp_dim(A22, 16 * GSL_DBL_EPSILON);
  1380.                     gsl_test(f, "  SV_decomp A(2x2)(%g, %g, %g, %g)", i1,i2,i3,i4);
  1381.                     s += f;
  1382.                   }
  1383.               }
  1384.           }
  1385.       }
  1386.   }
  1387.  
  1388.   return s;
  1389. }
  1390.  
  1391.  
  1392. int
  1393. test_cholesky_solve_dim(const gsl_matrix * m, const double * actual, double eps)
  1394. {
  1395.   int s = 0;
  1396.   size_t i, dim = m->size1;
  1397.  
  1398.   gsl_vector * rhs = gsl_vector_alloc(dim);
  1399.   gsl_matrix * u  = gsl_matrix_alloc(dim,dim);
  1400.   gsl_vector * x = gsl_vector_calloc(dim);
  1401.   gsl_matrix_memcpy(u,m);
  1402.   for(i=0; i<dim; i++) gsl_vector_set(rhs, i, i+1.0);
  1403.   s += gsl_linalg_cholesky_decomp(u);
  1404.   s += gsl_linalg_cholesky_solve(u, rhs, x);
  1405.   for(i=0; i<dim; i++) {
  1406.     int foo = check(gsl_vector_get(x, i), actual[i], eps);
  1407.     if(foo) {
  1408.       printf("%3d[%d]: %22.18g   %22.18g\n", dim, i, gsl_vector_get(x, i), actual[i]);
  1409.     }
  1410.     s += foo;
  1411.   }
  1412.   gsl_vector_free(x);
  1413.   gsl_matrix_free(u);
  1414.   gsl_vector_free(rhs);
  1415.  
  1416.   return s;
  1417. }
  1418.  
  1419. int test_cholesky_solve(void)
  1420. {
  1421.   int f;
  1422.   int s = 0;
  1423.  
  1424.   f = test_cholesky_solve_dim(hilb2, hilb2_solution, 2 * 8.0 * GSL_DBL_EPSILON);
  1425.   gsl_test(f, "  cholesky_solve hilbert(2)");
  1426.   s += f;
  1427.  
  1428.   f = test_cholesky_solve_dim(hilb3, hilb3_solution, 2 * 64.0 * GSL_DBL_EPSILON);
  1429.   gsl_test(f, "  cholesky_solve hilbert(3)");
  1430.   s += f;
  1431.  
  1432.   f = test_cholesky_solve_dim(hilb4, hilb4_solution, 2 * 1024.0 * GSL_DBL_EPSILON);
  1433.   gsl_test(f, "  cholesky_solve hilbert(4)");
  1434.   s += f;
  1435.  
  1436.   f = test_cholesky_solve_dim(hilb12, hilb12_solution, 0.5);
  1437.   gsl_test(f, "  cholesky_solve hilbert(12)");
  1438.   s += f;
  1439.  
  1440.   return s;
  1441. }
  1442.  
  1443.  
  1444. int
  1445. test_cholesky_decomp_dim(const gsl_matrix * m, double eps)
  1446. {
  1447.   int s = 0;
  1448.   size_t i,j, M = m->size1, N = m->size2;
  1449.  
  1450.   gsl_matrix * v  = gsl_matrix_alloc(M,N);
  1451.   gsl_matrix * a  = gsl_matrix_alloc(M,N);
  1452.   gsl_matrix * l  = gsl_matrix_alloc(M,N);
  1453.   gsl_matrix * lt  = gsl_matrix_alloc(N,N);
  1454.  
  1455.   gsl_matrix_memcpy(v,m);
  1456.  
  1457.   s += gsl_linalg_cholesky_decomp(v);
  1458.   
  1459.   /* Compute L LT */
  1460.   
  1461.   for (i = 0; i < N ; i++)
  1462.     {
  1463.       for (j = 0; j < N; j++)
  1464.         {
  1465.           double vij = gsl_matrix_get(v, i, j);
  1466.           gsl_matrix_set (l, i, j, i>=j ? vij : 0);
  1467.           gsl_matrix_set (lt, i, j, i<=j ? vij : 0);
  1468.         }
  1469.     }
  1470.             
  1471.   /* compute a = l lt */
  1472.   gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, l, lt, 0.0, a);
  1473.  
  1474.   for(i=0; i<M; i++) {
  1475.     for(j=0; j<N; j++) {
  1476.       double aij = gsl_matrix_get(a, i, j);
  1477.       double mij = gsl_matrix_get(m, i, j);
  1478.       int foo = check(aij, mij, eps);
  1479.       if(foo) {
  1480.         printf("(%3d,%3d)[%d,%d]: %22.18g   %22.18g\n", M, N, i,j, aij, mij);
  1481.       }
  1482.       s += foo;
  1483.     }
  1484.   }
  1485.  
  1486.   gsl_matrix_free(v);
  1487.   gsl_matrix_free(a);
  1488.   gsl_matrix_free(l);
  1489.   gsl_matrix_free(lt);
  1490.  
  1491.   return s;
  1492. }
  1493.  
  1494. int test_cholesky_decomp(void)
  1495. {
  1496.   int f;
  1497.   int s = 0;
  1498.  
  1499.   f = test_cholesky_decomp_dim(hilb2, 2 * 8.0 * GSL_DBL_EPSILON);
  1500.   gsl_test(f, "  cholesky_decomp hilbert(2)");
  1501.   s += f;
  1502.  
  1503.   f = test_cholesky_decomp_dim(hilb3, 2 * 64.0 * GSL_DBL_EPSILON);
  1504.   gsl_test(f, "  cholesky_decomp hilbert(3)");
  1505.   s += f;
  1506.  
  1507.   f = test_cholesky_decomp_dim(hilb4, 2 * 1024.0 * GSL_DBL_EPSILON);
  1508.   gsl_test(f, "  cholesky_decomp hilbert(4)");
  1509.   s += f;
  1510.  
  1511.   f = test_cholesky_decomp_dim(hilb12, 2 * 1024.0 * GSL_DBL_EPSILON);
  1512.   gsl_test(f, "  cholesky_decomp hilbert(12)");
  1513.   s += f;
  1514.  
  1515.   return s;
  1516. }
  1517.  
  1518.  
  1519. int
  1520. test_HH_solve_dim(const gsl_matrix * m, const double * actual, double eps)
  1521. {
  1522.   int s = 0;
  1523.   size_t i, dim = m->size1;
  1524.  
  1525.   gsl_permutation * perm = gsl_permutation_alloc(dim);
  1526.   gsl_matrix * hh  = gsl_matrix_alloc(dim,dim);
  1527.   gsl_vector * d = gsl_vector_alloc(dim);
  1528.   gsl_vector * x = gsl_vector_alloc(dim);
  1529.   gsl_matrix_memcpy(hh,m);
  1530.   for(i=0; i<dim; i++) gsl_vector_set(x, i, i+1.0);
  1531.   s += gsl_linalg_HH_svx(hh, x);
  1532.   for(i=0; i<dim; i++) {
  1533.     int foo = check(gsl_vector_get(x, i),actual[i],eps);
  1534.     if( foo) {
  1535.       printf("%3d[%d]: %22.18g   %22.18g\n", dim, i, gsl_vector_get(x, i), actual[i]);
  1536.     }
  1537.     s += foo;
  1538.   }
  1539.   gsl_vector_free(x);
  1540.   gsl_vector_free(d);
  1541.   gsl_matrix_free(hh);
  1542.   gsl_permutation_free(perm);
  1543.  
  1544.   return s;
  1545. }
  1546.  
  1547. int test_HH_solve(void)
  1548. {
  1549.   int f;
  1550.   int s = 0;
  1551.  
  1552.   f = test_HH_solve_dim(hilb2, hilb2_solution, 8.0 * GSL_DBL_EPSILON);
  1553.   gsl_test(f, "  HH_solve hilbert(2)");
  1554.   s += f;
  1555.  
  1556.   f = test_HH_solve_dim(hilb3, hilb3_solution, 128.0 * GSL_DBL_EPSILON);
  1557.   gsl_test(f, "  HH_solve hilbert(3)");
  1558.   s += f;
  1559.  
  1560.   f = test_HH_solve_dim(hilb4, hilb4_solution, 2.0 * 1024.0 * GSL_DBL_EPSILON);
  1561.   gsl_test(f, "  HH_solve hilbert(4)");
  1562.   s += f;
  1563.  
  1564.   f = test_HH_solve_dim(hilb12, hilb12_solution, 0.5);
  1565.   gsl_test(f, "  HH_solve hilbert(12)");
  1566.   s += f;
  1567.  
  1568.   f = test_HH_solve_dim(vander2, vander2_solution, 8.0 * GSL_DBL_EPSILON);
  1569.   gsl_test(f, "  HH_solve vander(2)");
  1570.   s += f;
  1571.  
  1572.   f = test_HH_solve_dim(vander3, vander3_solution, 64.0 * GSL_DBL_EPSILON);
  1573.   gsl_test(f, "  HH_solve vander(3)");
  1574.   s += f;
  1575.  
  1576.   f = test_HH_solve_dim(vander4, vander4_solution, 1024.0 * GSL_DBL_EPSILON);
  1577.   gsl_test(f, "  HH_solve vander(4)");
  1578.   s += f;
  1579.  
  1580.   f = test_HH_solve_dim(vander12, vander12_solution, 0.05);
  1581.   gsl_test(f, "  HH_solve vander(12)");
  1582.   s += f;
  1583.  
  1584.   return s;
  1585. }
  1586.  
  1587.  
  1588. int
  1589. test_TD_solve_dim(size_t dim, double d, double od, const double * actual, double eps)
  1590. {
  1591.   int s = 0;
  1592.   size_t i;
  1593.  
  1594.   gsl_vector * offdiag = gsl_vector_alloc(dim-1);
  1595.   gsl_vector * diag = gsl_vector_alloc(dim);
  1596.   gsl_vector * rhs = gsl_vector_alloc(dim);
  1597.   gsl_vector * x = gsl_vector_alloc(dim);
  1598.  
  1599.   for(i=0; i<dim; i++) {
  1600.     gsl_vector_set(diag, i, d);
  1601.     gsl_vector_set(rhs,  i, i + 1.0);
  1602.   }
  1603.   for(i=0; i<dim-1; i++) {
  1604.     gsl_vector_set(offdiag, i, od);
  1605.   }
  1606.  
  1607.   s += gsl_linalg_solve_symm_tridiag(diag, offdiag, rhs, x);
  1608.  
  1609.   for(i=0; i<dim; i++) {
  1610.     double si = gsl_vector_get(x, i);
  1611.     double ai = actual[i];
  1612.     int foo = check(si, ai, eps);
  1613.     if(foo) {
  1614.       printf("%3d[%d]: %22.18g   %22.18g\n", dim, i, gsl_vector_get(x, i), actual[i]);
  1615.     }
  1616.     s += foo;
  1617.   }
  1618.  
  1619.   gsl_vector_free(x);
  1620.   gsl_vector_free(rhs);
  1621.   gsl_vector_free(diag);
  1622.   gsl_vector_free(offdiag);
  1623.  
  1624.   return s;
  1625. }
  1626.  
  1627.  
  1628. int test_TD_solve(void)
  1629. {
  1630.   int f;
  1631.   int s = 0;
  1632.   double actual[16];
  1633.  
  1634.   actual[0] =  0.0;
  1635.   actual[1] =  2.0;
  1636.   f = test_TD_solve_dim(2, 1.0, 0.5, actual, 8.0 * GSL_DBL_EPSILON);
  1637.   gsl_test(f, "  solve_TD dim=2 A");
  1638.   s += f;
  1639.  
  1640.   actual[0] =  3.0/8.0;
  1641.   actual[1] =  15.0/8.0;
  1642.   f = test_TD_solve_dim(2, 1.0, 1.0/3.0, actual, 8.0 * GSL_DBL_EPSILON);
  1643.   gsl_test(f, "  solve_TD dim=2 B");
  1644.   s += f;
  1645.  
  1646.   actual[0] =  5.0/8.0;
  1647.   actual[1] =  9.0/8.0;
  1648.   actual[2] =  2.0;
  1649.   actual[3] =  15.0/8.0;
  1650.   actual[4] =  35.0/8.0;
  1651.   f = test_TD_solve_dim(5, 1.0, 1.0/3.0, actual, 8.0 * GSL_DBL_EPSILON);
  1652.   gsl_test(f, "  solve_TD dim=5");
  1653.   s += f;
  1654.  
  1655.   return s;
  1656. }
  1657.  
  1658. int
  1659. test_bidiag_decomp_dim(const gsl_matrix * m, double eps)
  1660. {
  1661.   int s = 0;
  1662.   size_t i,j,k,r, M = m->size1, N = m->size2;
  1663.  
  1664.   gsl_matrix * A  = gsl_matrix_alloc(M,N);
  1665.   gsl_matrix * a  = gsl_matrix_alloc(M,N);
  1666.   gsl_matrix * b  = gsl_matrix_alloc(N,N);
  1667.  
  1668.   gsl_matrix * u  = gsl_matrix_alloc(M,N);
  1669.   gsl_matrix * v  = gsl_matrix_alloc(N,N);
  1670.  
  1671.   gsl_vector * tau1  = gsl_vector_alloc(N);
  1672.   gsl_vector * tau2  = gsl_vector_alloc(N-1);
  1673.   gsl_vector * d  = gsl_vector_alloc(N);
  1674.   gsl_vector * sd  = gsl_vector_alloc(N-1);
  1675.  
  1676.   gsl_matrix_memcpy(A,m);
  1677.  
  1678.   s += gsl_linalg_bidiag_decomp(A, tau1, tau2);
  1679.   s += gsl_linalg_bidiag_unpack(A, tau1, u, tau2, v, d, sd);
  1680.  
  1681.   gsl_matrix_set_zero(b);
  1682.   for (i = 0; i < N; i++) gsl_matrix_set(b, i,i, gsl_vector_get(d,i));
  1683.   for (i = 0; i < N-1; i++) gsl_matrix_set(b, i,i+1, gsl_vector_get(sd,i));
  1684.   
  1685.   /* Compute A = U B V^T */
  1686.   
  1687.   for (i = 0; i < M ; i++)
  1688.     {
  1689.       for (j = 0; j < N; j++)
  1690.         {
  1691.           double sum = 0;
  1692.  
  1693.           for (k = 0; k < N; k++)
  1694.             {
  1695.               for (r = 0; r < N; r++)
  1696.                 {
  1697.                   sum += gsl_matrix_get(u, i, k) * gsl_matrix_get (b, k, r)
  1698.                     * gsl_matrix_get(v, j, r);
  1699.                 }
  1700.             }
  1701.           gsl_matrix_set (a, i, j, sum);
  1702.         }
  1703.     }
  1704.  
  1705.   for(i=0; i<M; i++) {
  1706.     for(j=0; j<N; j++) {
  1707.       double aij = gsl_matrix_get(a, i, j);
  1708.       double mij = gsl_matrix_get(m, i, j);
  1709.       int foo = check(aij, mij, eps);
  1710.       if(foo) {
  1711.         printf("(%3d,%3d)[%d,%d]: %22.18g   %22.18g\n", M, N, i,j, aij, mij);
  1712.       }
  1713.       s += foo;
  1714.     }
  1715.   }
  1716.  
  1717.   gsl_matrix_free(A);
  1718.   gsl_matrix_free(a);
  1719.   gsl_matrix_free(u);
  1720.   gsl_matrix_free(v);
  1721.   gsl_matrix_free(b);
  1722.   gsl_vector_free(tau1);
  1723.   gsl_vector_free(tau2);
  1724.   gsl_vector_free(d);
  1725.   gsl_vector_free(sd);
  1726.  
  1727.   return s;
  1728. }
  1729.  
  1730. int test_bidiag_decomp(void)
  1731. {
  1732.   int f;
  1733.   int s = 0;
  1734.  
  1735.   f = test_bidiag_decomp_dim(m53, 2 * 64.0 * GSL_DBL_EPSILON);
  1736.   gsl_test(f, "  bidiag_decomp m(5,3)");
  1737.   s += f;
  1738.  
  1739.   f = test_bidiag_decomp_dim(m97, 2 * 64.0 * GSL_DBL_EPSILON);
  1740.   gsl_test(f, "  bidiag_decomp m(9,7)");
  1741.   s += f;
  1742.  
  1743.   f = test_bidiag_decomp_dim(hilb2, 2 * 8.0 * GSL_DBL_EPSILON);
  1744.   gsl_test(f, "  bidiag_decomp hilbert(2)");
  1745.   s += f;
  1746.  
  1747.   f = test_bidiag_decomp_dim(hilb3, 2 * 64.0 * GSL_DBL_EPSILON);
  1748.   gsl_test(f, "  bidiag_decomp hilbert(3)");
  1749.   s += f;
  1750.  
  1751.   f = test_bidiag_decomp_dim(hilb4, 2 * 1024.0 * GSL_DBL_EPSILON);
  1752.   gsl_test(f, "  bidiag_decomp hilbert(4)");
  1753.   s += f;
  1754.  
  1755.   f = test_bidiag_decomp_dim(hilb12, 2 * 1024.0 * GSL_DBL_EPSILON);
  1756.   gsl_test(f, "  bidiag_decomp hilbert(12)");
  1757.   s += f;
  1758.  
  1759.   return s;
  1760. }
  1761.  
  1762.  
  1763. int main(void)
  1764. {
  1765.   gsl_ieee_env_setup ();
  1766.  
  1767.   m35 = create_general_matrix(3,5);
  1768.   m53 = create_general_matrix(5,3);
  1769.   m97 = create_general_matrix(9,7);
  1770.  
  1771.   s35 = create_singular_matrix(3,5);
  1772.   s53 = create_singular_matrix(5,3);
  1773.  
  1774.   hilb2 = create_hilbert_matrix(2);
  1775.   hilb3 = create_hilbert_matrix(3);
  1776.   hilb4 = create_hilbert_matrix(4);
  1777.   hilb12 = create_hilbert_matrix(12);
  1778.  
  1779.   vander2 = create_vandermonde_matrix(2);
  1780.   vander3 = create_vandermonde_matrix(3);
  1781.   vander4 = create_vandermonde_matrix(4);
  1782.   vander12 = create_vandermonde_matrix(12);
  1783.  
  1784.   moler10 = create_moler_matrix(10);
  1785.  
  1786.   c7 = create_complex_matrix(7);
  1787.  
  1788.   row3 = create_row_matrix(3,3);
  1789.   row5 = create_row_matrix(5,5);
  1790.   row12 = create_row_matrix(12,12);
  1791.  
  1792.   A22 = create_2x2_matrix (0.0, 0.0, 0.0, 0.0);
  1793.  
  1794.   /* Matmult now obsolete */
  1795. #ifdef MATMULT
  1796.   gsl_test(test_matmult(),        "Matrix Multiply"); 
  1797.   gsl_test(test_matmult_mod(),    "Matrix Multiply with Modification"); 
  1798. #endif
  1799.   gsl_test(test_bidiag_decomp(),  "Bidiagonal Decomposition");
  1800.   gsl_test(test_LU_solve(),       "LU Decomposition and Solve");
  1801.   gsl_test(test_LUc_solve(),      "Complex LU Decomposition and Solve");
  1802.   gsl_test(test_QR_decomp(),      "QR Decomposition");
  1803.   gsl_test(test_QR_solve(),       "QR Solve");
  1804.   gsl_test(test_QR_lssolve(),     "QR LS Solve");
  1805.   gsl_test(test_QR_update(),      "QR Rank-1 Update");
  1806.   gsl_test(test_QRPT_decomp(),    "QRPT Decomposition");
  1807.   gsl_test(test_QRPT_solve(),     "QRPT Solve");
  1808.   gsl_test(test_SV_decomp(),      "Singular Value Decomposition");
  1809.   gsl_test(test_SV_solve(),       "SVD Solve");
  1810.   gsl_test(test_cholesky_decomp(),"Cholesky Decomposition");
  1811.   gsl_test(test_cholesky_solve(), "Cholesky Solve");
  1812.   gsl_test(test_HH_solve(),       "Householder solve");
  1813.   gsl_test(test_TD_solve(),       "Tridiagonal solve");
  1814.  
  1815.   gsl_matrix_free(hilb2);
  1816.   gsl_matrix_free(hilb3);
  1817.   gsl_matrix_free(hilb4);
  1818.   gsl_matrix_free(hilb12);
  1819.  
  1820.   gsl_matrix_free(vander2);
  1821.   gsl_matrix_free(vander3);
  1822.   gsl_matrix_free(vander4);
  1823.   gsl_matrix_free(vander12);
  1824.  
  1825.   gsl_matrix_complex_free(c7);
  1826.  
  1827.   exit (gsl_test_summary());
  1828. }
  1829.